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We refit the NRL tight binding parameterization for Aluminium by Mehl et al [Phys. Rev. B, 6f, 
4894 (2000)], to a database generated via full potential Linearized Augmented Plane Wave (LAPW) 
Density Functional Theory (DFT) calculations. This is performed using a global optimization 
algorithm paying particular attention to reproducing the correct order of the angular symmetries 
of the tight binding fee and bec bandstructure. The resulting parameterization is found to better 
predict the hep phase and both the stable and unstable planar stacking fault defect energies. 
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The empirical tight binding technique can be seen as 
a compromise between the very accurate but computa- 
tionally expensive ah initio Density Functional Theory 
(DFT) methods, and the fast but less accurate empirical 
potential methods such as the Embedded Atom Method 
(EAM). Although great reductions in computational cost 
has been achieved for ah initio methods by employing 
a minimal basis set and efficient algorithms that scale 
linearly with atom number 1 , fully self-consistent calcula- 
tions are still relatively slow for "large" scale Molecular 
Dynamics (MD) simulations. For empirical tight bind- 
ing methods, a range of linear scaling techniques have 
also been developed^ and the development of accurate 
tight binding parameterizations is of utmost importance, 
if this method is to remain competitive. Furthermore, the 
strength of the tight binding method, contrary to many 
EAM parameterizations, is its ability to accurately pre- 
dict a range of structural properties not included in the 
fitting database. 

The NRL tight binding formalism is a total energy 
tight binding method based on the fitting of bandstruc- 
ture eigenvalues and total energies from ab initio calcu- 
lations. One of the advantages of this method is that one 
is able to fit to the bandstructure and total energy of a 
given structure simultaneously be redefining the energy 
zero of the sum of eigenvalues. From DFT, the total en- 
ergy can be expressed as the sum of the single electron 
eigenvalues plus an ionic-like term 



rection term" 



E(n(v)) 



/(p-e<)ei + F(n(r)) 



(1) 



where f{x) is the Fermi function, ej the electronic eigen- 
values, /i the chemical potential, and F(n(r)) the so- 
called "correction term" 3,4 . Following Mehl et al this 
can be reformulated by shifting the energy zero of the 
bandstructure term by a value exactly canceling the "cor- 



£(n(r))=£/(y-e^ 
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This shift is taken into account by making the onsite 
terms environmentally dependent. Using a nonorthogo- 
nal two-center Slater-Koster formulation^, the tight bind- 
ing model is then parameterized by fitting the bandstruc- 
ture and total energies as a function of lattice constant 
and crystal structure of the tight binding model to a 
database produced by ab inito calculations. The philoso- 
phy of the NRL formalism is to include only a handful of 
high symmetry structures in the fitting database. This 
method has proven to be successful for a wide range of 
elements. Further details can be found in Refsi&ii. 

We have refitted a tight binding parameterization by 
Mehl et a& for Aluminium using an s, p and d or- 
bital basis, to the full potential LAPW DFT calcula- 
tions of WIEN2k^. As an alternative to the conventional 
least squares procedure we have used Adaptive Simu- 
lated Annealing (ASA), a stochastic global optimization 
algorithrrAiS. Assuming a probability density of states 
<7t(x) for the parameters space {xi}, and a cost function 
defined by the "energy" E(x). At any step, k, the "prob- 
ability" for accepting the new cost function E^+i is given 
by: 



h(AE) 



exp(-AE/T) for AE > 0, 



1 



for AE < 0. 



(3) 



where the parameter T = T{k) defines the "annealing 
schedule". For such a system the principle of detailed 
balance holds, which means that, in theory, all states of 
the system will be sampled by the algorithm and a global 
minimum can be found. A weakness of the simulated 
annealing approach is that it needs to sample a large 
part of parameter space and can therefore be inherently 
slow. This is partly remedied in the ASA algorithm by 
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assuming a probability density for the parameters with 
"fat" logarithmic tails 

9T ( X ) s RAM = n 2(N+ ^ ln(1+7V) ( 4 ) 

where x = {xi G [— 1, 1]; i = 1, N} is the normalized pa- 
rameter set. This allows for a very fast annealing sched- 
ule 

Ti(fc) =T to exp-c i fe 1 / JV (5) 

in contrast to the conventional linearly decreasing anneal- 
ing schedule. Tio and Cj are free parameters, which allows 
the algorithm to adapt to a changing "environment" for 
the cost functions as the search proceeds. Further details 
can be found in ref. @- 

In the fitting database we included band-structure data 
as a function of lattice constant for the the fee and bec 
crystal structures using irreducible k-point meshes of 146 
for fee and 91 points for bee. In order to ensure the 
correct angular character for each eigenvalue and fur- 
ther constrain the Slater-Koster parameters we block- 
diagonalized the Hamiltonian for a large number of high 
symmetry points and directions: the T, X, L, W points 
and A, A, E, S, Z directions for the fee structure, and 
the r, H, P, N points and A, A, E, F, D, G directions 
for the bec structure. 

For the DFT calculations we used the input parame- 
ters: R m t = 2.5 Bohr, R m tK max = 7 and G max = 14 
Ry 1 / 2 . Here, R m t is the muffin tin radius, K max is the 
plane wave cut-off and G max is the maximum Fourier 
component of the electron density. For the exchange- 
correlation potential we used the Generalized Gradient 
Approximation (GGA) of Perdew et. alii. 

Fig-ffldisplays the TB bandstructure for fee at a = 7.65 
Bohr compared to FLAPW ab inito calculations and the 
TB parameterization of Mehl et. al. In contrast to the 
previous parametrization, the ordering with respect to 
angular character of all bands below 1.25 Ry is now cor- 
rect. This is important, since otherwise the relative con- 
tribution to the bandstructure from each atomic orbital 
in the LCAO-expansion would break the basic symmetry- 
requirements of the crystal potential^. 

Fig. [3 displays the energies of state for the fee, hep and 
bec structures as obtained from a second order Birch fit 
to the data points derived from the present TB model. 
Also shown are FLAPW ab inito data for fee and bec. 
Tab. Q] compares the elastic constants from the TB pa- 
rameterization of the present work with the parameter- 
ization of Mehl et aim, ab initio DFT calculations using 
WIEN2k and experimental values. As can be seen, the 
TB values are in excellent agreement with the ab inito 
and experimental values, where only the bulk modulus is 
overestimated. This later discrepancy is probably due to 
more weight being given (in the present fitting process) 
to bandstructure data than the equation of state total 
energy data. The hep data is in effect a prediction of the 
tight binding parameterization since only the fee and bec 
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Band Structure of Al from TB (New) 
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b) 

Al Band Structure from WIEN2k 
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Band Structure of Al from TB (MP) 
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FIG. 1: Band structure for fee Al at a = 7.65 Bohr for a) the 
new TB parameterization, b) WIEN2k ab initio calculations, 
and c) the TB parameterization of Mehl et. al. 



structures where fitted. At the obtained equilibrium lat- 
tice constants, our fit predicts a fec-hep energy difference 
of 2.1 mRyd per atom, which agrees well with that of 
the converged (with respect to k-point mesh density us- 
ing the previously stated ab initio parameters) WIEN2k 
value of 2.0 mRyd. The TB parametrization of Mehl et. 
al. returns a value of 1.8 mRyd. 

An important fee material property for any interatomic 
model to reproduce is the intrinsic stacking fault surface 
energy density since it plays a crucial role in the disas- 
sociation width of a full dislocation into its constituent 
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FIG. 2: Energies of state for the bec, fee, and hep phases. 
Note that only information from bec and fee structures was 
included in the fitting database. 




FIG. 3: Bain Path at a — 7.65 Bohr comparing the new 
parameterization with that of Mehl et. al (MP). 



TABLE I: Elastic constants for fee in units of GPa calculated 
at a = 7.65 Bohr. Experimental results are taken from ref. 



Property 


TB 


TB MP 


DFT 


Exp. 


a 


7.64 


7.58 


7.65 


7.65 


B 


92 


79 


78 


79 


Cll — Cl2 


48 


59 


50 


46 


C44 


27 


23 


27 


28 



Shockley partials and is believed to play a role in the is- 
sue of why in some cases only partial dislocations are ob- 
served in atomistic modeling of the mechanical properties 
of nanocrystalline materials 14,15 . An unrelaxed stack- 
ing fault is modeled by stacking layers of closed packed 
planes in the [111] direction using a particular stacking 
sequence. We choose the primitive vectors of the unit 
cell to be: 

1 1 
ai =-ay + -az 

& 2 =2 ax + 2° z ^ 

as =(4 + |)ax + (4 + |)ay - (4 - |)az 
6 6 3 

where a is the lattice constant, q is a variable controlling 
the slip in the [112] direction. This is equivalent to a su- 
percell of 12 close-packed planes, with a variable slip at 
the 12th layer due to the periodic boundary conditions. 
For a q = we get a perfect fee crystal represented by 
the stacking sequence ABC ABC, and for q = 1 we get 
the full stacking fault ABC\BC, where the stacking se- 
quence at the interface corresponds to two 111 planes 
in which the atoms are hep coordinated. A mapping of 
the structural energy for all values of q from to 1 is 
called a Bain Path (or also a generalized stacking fault 
curve) where an intrinsic staking fault occurs at q = 1. 
Fig. shows the stacking fault energy as a function of 
the slip variable, q relative to the fee energy at q = 0. 
The maximum of the curve defines the unstable stacking 



fault configuration, which in fig. [31 is near q = 0.6. 

Table [n] presents the calculated values for the intrin- 
sic stacking fault energy, ji s , unstable stacking fault en- 
ergy, 7„ s , and the ratio Ji s /jus, compared to the TB pa- 
rameterization of Mehl et. ali&. In addition we include 
published results from DFT and EAM calculations, and 
experimental values. The experimentally derived values 
for 7i s vary from a low value of 110 mJ/m 2 to a high 
value of 280 mJ/m 2 . However, ab initio DFT calcula- 
tions have consistently shown results in the range 161-166 
mJ/m 2 . The EAM potential of Voter and Chen shows 
values for the stacking fault energies that are seriously 
underestimated. The important point here is that these 
values were not included in the fitting database for this 
potential^. Similar results have been produced by the 
EAM potential of Ercolessi and Adams, which has been 
used in large scale Molecular Dynamics (MD) simula- 
tions of extended dislocations 14 . In contrast to this, the 
recently published EAM potential of Mishin and Farkas, 
show values for the intrinsic and unstable stacking fault 
energies which are more in line with the ab initio cal- 
culations. The reason for this is that a value for the 
intrinsic stacking fault energy was included in the fit- 
ting database^. The TB parameterization of the present 
work shows a marked improvement over the older param- 
eterization of Mehl et al Especially the intrinsic stacking 
fault energy, ji s , with a value of 146 mJ/m 2 which is in 
much better agreement with ab inito calculations than a 
value of 97 mJ/m 2 of the original fit. Another important 
improvement is the value of the ratio jus /lis, which is 
a measure of the relative cost of dislocation nucleation 
with respect to the intrinsic stacking fault energy. For 
the new TB parameterization this value is 1.3 which is 
in good agreement with the value of 1.4 from ab inito 
calculations. 

The improved prediction of the stacking fault energies 
is consistent with the corresponding improvement of the 
hep cohesive energy and lattice constant of the present 
model. At the nearest neighbor bond length, the local fee 
and hep environments differ only with respect to their di- 
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TABLE II: Unstable and intrinsic stacking fault energies for Aluminium in mj/m 2 



Property 


TB 

Present work 


TB 
Mehl et al a 


ab initio 
DFT 6 


EAM 
Mishin and Farkas c 


EAM 
Voter and Chen d 


Exp. c 


7« s 


183 


177 


224 


167 


93 






146 


97 


166 


146 


76 


110-216 




1.3 


1.8 


1.4 


1.2 


1.2 





"Usin g th e parameters of ref. Ilfi 
6 Ref. p] 
c Ref. J| 
d Ref. Ilj 

e High value from ref. l20t low value ref. [2l] 



hedral angle arrangements. In terms of energy and lattice 
constant this difference should be sensitive to the angular 
character of, particularly, the n and S bonds of the p and 
d states, and it is precisely such bonding elements that 
the current parameter-set more accurately reproduces at 
and near the Fermi level. 

In conclusion, we have refitted the tight binding pa- 
rameters of Mehl et al for Aluminium to ab initio 
FLAPW eigenvalues. We have used a stochastic global 
optimization algorithm and an extended database of high 
symmetry eigenvalues. As a result we obtain a parame- 



terization with an improved representation of the angular 
symmetries for the Slater Koster parameters. As a con- 
sequence we obtain improved predictions for the stacking 
fault energies without including this information explic- 
itly in the fitting database. 
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